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Abstract 

We study a theoretical model of mud cracks, that is, the fracture patterns 

resulting from the contraction with drying in a thin layer of a mixture of 

granules and water. In this model, we consider the slip on the bottom of this 

layer and the relaxation of the elastic field that represents deformation of the 

layer. Analysis of the one-dimensional model gives results for the crack size 

that are consistent with experiments. We propose an analytical method of 

estimation for the growth velocity of a simple straight crack to explain the 

very slow propagation observed in actual experiments. Numerical simulations 

reveal the dependence of qualitative nature of the formation of crack patterns 

on material properties. 

46.35. +z,46.50.+a,47.54.+r,62.20.Mk 



Typeset using REVTeX 



1 kit sune@minnie . disney. phys . nar a- wu. ac . j p 



1 



I. INTRODUCTION 



Many kinds of mixtures of granules and water, such as clay, contract upon desiccation 
and form cracks. These fracture patterns are familiar to us as ordinary mud cracks. However, 
the fundamental questions about these phenomena have not yet been answered theoretically. 
The problems which need to be addressed include determining the condition under which 
fragmentation occurs, the dynamics displayed by cracks, and the patterns which grow. 

In simple and traditional experiments on mud cracks, a thin layer of a mixture in a 
rigid container with a horizontal bottom is prepared left to dry at room temperature |T] -Q. 
Typically, clay, soil, flour, granules of magnesium carbonate and alumina are used. In 
almost all cases, cracks extend from the surface to the bottom of the layer and propagate 
horizontally along a line, forming a quasi-two-dimensional structure. Typically we observe 
a tiling pattern composed of rectangular cells in which cracks mainly join in a T-shape. 
Groisman and Kaplan carried out more detail experiments with coffee powder and reported 
1) that the size of a crack cell after full drying is nearly proportional to the thickness of 
the layer and larger in the case of a "slippery" bottom, 2) that the velocity of a moving 
crack is almost independent of time for a given crack and very slow on the order of several 
millimeters per minute, but that it differs widely from one crack to another, and 3) that as 
the layer becomes thin, there is a transition to patterns which contain many Y-shape joints 
and unclosed cells owing to the arrest of cracks [[§. 

Another experimental setup was used by Allain and Limat |J . This setup produces cracks 
that grow directionally by causing evaporation to proceed from one side of the container. 
Sasa and Komatsu have proposed a theoretical model for such systems [0. 

Fragmentation of coating or painting also arises from desiccation, This has been studied 
theoretically by some people pHlOfl. From the viewpoint that fractures are caused by slow 
contraction, these problems can be thought of as belonging to the same category as thermal 
cracks in glasses [pH lOU and the formation of joints in rocks brought by cooling ||T7| , |T8 



In addition, we note that mixtures of granular matter and fluid have properties that vary 
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greatly from that of complete elastic materials, in particular, dissipation and viscoelasticity. 
The propagation of cracks in such media has been investigated recently using developments 
in nonlinear physics [p~9| |28[] . 



In this paper, we undertake a theoretical investigation of the experiments described 
above. We treat such system as consisting of fractures arising from quasi-static and uniform 
contraction in thin layers of linear elastic material. 

In Sec. [ll|, we propose a one-dimensional model. Our model takes into account the slip 
displacement on the bottom of a container, because most of the experiments can not be 
assumed to obey a fixed boundary condition. We can investigate the development of the 
size of a crack cell by applying a fragmentation condition to the balanced states of the elastic 
field. 



In Sec. Ill], we report the analytical results of our one-dimensional model. Here we 
consider both the critical stress condition and the Griffith criterion as the fragmentation 
condition. We consider these two alternative criteria because the nature of the breaking 
condition in mixtures of granules and water is not clear. The critical stress condition predicts 
that the final size of a crack cell is proportional to the thickness of the layer and that, in 
the case of a slippery bottom, it becomes much larger than the thickness. These predictions 
seem to be consistent with the experimental results. In contrast, we find that the Griffith 
criterion predicts a different relation between the final size of a cell and the thickness. 

In Sec. |V|, we extend the model to two dimensions and investigate the time development 
of a crack. In order to describe the relaxation process of the elastic field, we use the Kelvin 
model while taking into account the effect of the bottom of the container. We assume the 
stress, excluding dissipative force, to be constant in the front of a propagating crack tip and 
evaluate the velocity of a simple straight crack tip analytically. Our results indicate that 
cracks advance at very slow speed in comparison with the sound velocity. 

In Sec. [V], we report on the numerical simulations of our model that reproduce fracture 
patterns similar to those in real experiments. The growth of the patterns exhibits qualitative 



differences depending on the elastic constants and the relaxation time. As the relaxation 
time becomes smaller, in particular, we observe the growth of fingering patterns with tip 
splitting rather than side branching of cracks. 

Finally, we conclude the paper with a summary of the results and a discussion of the 
open problems in Sec. |V|. 

II. MODELING OF FRACTURE CAUSED BY SLOW SHRINKING 

We analyze the formation of cracks induced by desiccation in terms of the following four 
processes. 

1. The water in a mixture evaporates from the surface of a layer. 

2. Each part of the mixture shrinks upon desiccation. 

3. Stress increases in the material because contraction is hindered near the bottom of a 
container. 

4. Fracture arises under some fragmentation condition. 

In this section, we examine each process individually and construct a one-dimensional 
model, where we introduce simple assumptions regarding the unclear properties of granu- 
lar materials. Some similar models have been proposed previously f|[|]]. One-dimensional 
models assume that cracks are formed one at a time, each propagating along a line and 
thereby dividing the system into two pieces separated by a boundary with one-dimensional 
structure. Using this assumption, we can ignore the propagation of cracks and consider the 
development of patterns by using only the condition of separation. 

1. From a microscopic viewpoint, water either exists in the inside of the particles of gran- 
ular materials or acts to create bonds between the particles. Here, we can introduce 
the water content averaged over a much larger area than that of a single particle and 
measure the degree of drying. When the thickness of a layer H is sufficiently thin and 
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the characteristic time of desiccation Td is very large, the water content in the layer 
can be considered uniform. Assuming that water transfers diffusively in a layer, the 
sufficient condition here is that H 2 /Td is much smaller than the diffusion constant. 
Therefore, we restrict our consideration to the case of the uniform water distribution 
and exclude the process of water transfer from the model. 

2. The main cause of contraction is the shrinking of particles in the mixture arising 
from desiccation. The water content is considered to determine the shrinking rate 
in the case of uniform contraction in which all the boundaries of the mixture are 
stress-free. We refer to this shrinking rate as "free shrinking rate" in the following 
discussions and this concept is used in place of the concept of the water This makes 
clear the relation between the present problems and those involving fractures induced 



by other causes, such as temperature gradient with slow contraction. We note, 

however, that it is more difficult to measure the free shrinking rate than the water 
content experimentally and it is thus necessary to know their relation to compare our 
theory with experiments on the time development of patterns. 

The contraction force is thought to arise from the water bonds among particles. We 
estimate the Reynolds number R e to consider the behavior of the water in a bond. The 
diameter of a particle R is generally about 0.1mm, and the kinematic viscosity of water 
v is about 1mm 2 js. Although the propagation of a crack causes the displacement of 
surrounding particles with opening the crack surfaces, the velocity of the displacement 
is smaller than the crack speed itself, except in the microscopic region at the crack tip. 
Therefore we estimate the typical velocity of water V in the bulk of a mixture to be 
smaller than the crack speed. The crack speed has been measured as about 0.1 mm/s 
in experiments and it is, of course, considerably faster than the shrinking speed of 
the horizontal boundary with desiccation, which is typically about 10mm/ day. Thus 
the Reynolds number R e = RV/v is estimated to be smaller than 1/100. We expect 
that the water among the particles behaves like a viscous fluid and that the mixture 
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displays strong dissipation. 

If a material displays strong dissipation and shrinks quasi-statically, the elastic field 
is balanced to minimize the free energy, except during the time when cracks are prop- 
agating. We deal with balanced states for the present and return to the problem of 
relaxation in order to treat the development of cracks in Sec. [TV]. Because mixtures of 
granules and fluids have many unclear properties with respect to elasticity, we idealize 
them as linear elastic materials. When the free volume-shrinking rate C v is uniform, 
it is well known that the free energy density of a uniform and isotropic linear elastic 



material is given in terms of the stress tensor Uij in the form |30 



1 / 1 \ 2 

e v = -k v (uu + C v ) 2 + fi \u ik - -uii5 ik j , (2.1) 

where n v and \i are the elastic constants and repeated indices indicate summation. 
The stress tensor is expressed as 



3e 

cr v ij = -t^- = (A«n + K v C v )5ij + 2/j.Uij, (2.2) 



where k v = A+2/z/3. As a result, the balanced equation of the elastic field da v ij/ dxj = 
does not include the shrinking rate C v for linear elastic materials with uniform 
contraction and it is the same as in the case of the elastic materials without shrinking. 
The effect of contraction appears only through the boundary conditions. 

After the formation of cracks divides the system into cells, each cell is independent 
of the others, because the vertical surfaces of cracks become stress-free boundaries. 
Without considering the boundary conditions on the lateral sides of the container, 
we can simplify the problem by starting with the initial condition that the system 
has stress-free boundaries on its lateral sides. In contrast to the lateral and the top 
surfaces, the bottom of a layer is not a stress-free boundary. The difference among 
the boundary conditions produces strain with contraction and then stress. This is the 
cause of fracture. 



We observe the slip of layers along the bottom in most experiments. Thus we introduce 
slip displacement with a frictional force into the model. Because the frictional force is 
caused by the water between the bottom of a layer and the container, it is considered 
to remain finite even in the limit of vanishing thickness of a layer if 0. In order to 
understand the effects of friction for crack patterns, we simplify the maximum frictional 
force per unit area of the surface to be constant and independent of H without making 
a distinction between static and kinetic frictional effects. 

4. Cracks propagate very slowly in a mixture of granules and water. Because this propa- 
gation resembles quasi-static growth of cracks, the first candidate of the fragmentation 
condition is the Griffith criterion applied to the free energy of the entire system. 

However, we note that ordinary brittle materials break instantaneously, not quasi- 
statically, in the situation that the stress increases without fixing the deformation 
of the system. The situation is similar to that in a shrinking mixture. We need to 
consider the possibility that cracks in a mixture propagate slowly owing to dissipation. 
Therefore we consider two typical fragmentation conditions, the critical stress condition 
and the Griffith criterion, in the first and the last halves of Sec.^TJ respectively. 

In the context of the critical stress condition, the fragmentation condition is that 
the maximum principal stress exceeds a material constant at breaking. This has also 
been used in many numerical models because of the technical advantage of the local 
condition. Our model introduces the critical value for the energy density as an equiv- 
alent condition. We note that the energy of the system before the fragmentation is 
higher than after the fragmentation because the critical value is a material constant 
independent of the system size. 

In contrast, using the Griffith criterion as an alternative fragmentation condition stip- 
ulates that the energy changes neither before nor after breaking. This condition is 
used by Sasa and Komatsu in their theory 
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With the above considerations, we construct a one-dimensional model, following the lead 
of Sasa and Komatsu [^]. As we show in Fig. |1], we consider a chain of springs by a distance 
a as the discrete model of the thin layer of an elastic material, where we number the nodes 
i = —N, —N + 1, N — 1, N for a system of half-size L := aN. In order to represent the 
vertical direction of a sufficiently thin layer, we introduce vertical springs with length H 
which connect each node to an element on the bottom. The vector (ui,Vi) represents the 
horizontal and vertical displacements of the ith node, and w% is the horizontal displacement 
of the element on the bottom connected with the node. The shrinking of a material is 
modeled by decreasing the natural length of the springs. We assume an isotropic material 
with a linear free shrinking rate s, making the natural lengths of both the horizontal and 
vertical springs to be 1 — s times the initial lengths, i.e. (1 — s)a and (1 — s)H. 

If the vertical springs are simple ordinary springs, linear response is lost under shearing 
strain. We therefore add non-simple springs to the vertical direction which produce a hor- 
izontal force in the case that Ui ^ Wi to represent a linear elastic material. This type of 
spring is used in the model of Hornig et al. |J . The energy of the system is described by 

j AT-l 

E = - K^Ui+i - ui + sa) 2 

^ i=~N 

1 N 

+ - £ [K 2 (u t - Wl ) 2 + K' 2 (vi + sH)% (2.3) 

1 i=—N 

where Ki,K 2 and K' 2 are the spring constants. Because f« is included in the last term 
independently of both it, and u>,, it follows that Vi = —sH in the balanced states, and then 
this term vanishes. This indicates that the model neglects the horizontal stress arising from 
vertical contraction. Thus all we need to do is minimize the energy ( |2.3j ) without the last 
term in order to find the horizontal displacements U{ and Wi. In the continuous limit, a — > 0, 
the above energy should be described using an independent energy density for H, as in the 
case of ( |2.1| ) for a linear elastic material. We scale the space length by the thickness of 
a layer H and introduce the space coordinate x := ai/H. Through the transformation to 
non-dimensional variables L — > LH, u% — > Hu(x), Wi — > Hw(x) and E — > H 2 E, the energy 
( ^73|) becomes 
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E = f dx{ ei (x) + e 2 (x)}, (2.4a) 
e 1 (x) = -k 1 (u x + s) 2 (2.4b) 
and e 2 (x) = ~k 2 (u — w) 2 , (2.4c) 

and we know that both k\ := Jfi^i and k 2 '■= —K 2 are the independent constants of H. 

We introduce the maximum frictional force F s for the slip of the elements on the bottom, 
as explained above. The vertical spring pulls the zth element along the bottom with the 
force Fi = K 2 {u{ —Wi). Each element on the bottom remains stationary if |i^| < F s , and, if 
not, it slips to a position at which = F s is satisfied. The slip condition is expressed by 
the energy density of a vertical spring through the following rule in the previous continuous 
a — > limit: 

1 s 

e 2 (x) > -hsl =^> w(x) = u(x) ± — . (2.5) 

2 q 

Here the choice of the sign depends on the direction of the force. The constants s s and q 
are defined by the equations 

5*^ s 2& and « s /i- (2 - 6) 

We note the constant q is order 1, because its square represents a ratio of certain elastic 
constants which are the same order in ordinary materials. 

Neglecting the short periods during which the system experiences cracking and slip, (|2.4|) 
and ( j2.5|) constitute the closed form of our one- dimensional model with the fragmentation 
condition given in the next section. 

III. ANALYSIS OF THE ONE-DIMENSIONAL MODEL 

We here report analytical results of our one-dimensional model for the typical two frag- 
mentation conditions, i.e, the critical stress condition and the Griffith condition. 



9 



A. The Critical Stress Condition 



The critical stress condition demands that the maximum principal stress exceeds a ma- 
terial constant at the fragmentation. 

In the case of this condition, we can generally demonstrate that it is difficult to treat 
the bottom surface as a fixed boundary for a uniform and isotropic elastic material. We 
first explain it before the analysis of the one-dimensional model. Let us think of the layer 
of a linear elastic material contracting with a fixed boundary condition on the bottom. It is 
shrinking more near the top surface, and the cross section assumes the form of a trapezoid 
as we show schematically in Fig. ^|. We compare the stress at the following three points: 
(A) the horizontal center of the cell near the bottom; (B) the lateral point near the bottom; 
and (C) the horizontal center above the bottom. The horizontal tensions at A and B are the 
same because of the fixed boundary condition on the bottom. Although the stress at C is as 
horizontal as at A, the strength is weaker. B is also pulled in the direction along the lateral 
surface due to the deformation. Using A, B, and C to represent the respective strengths 
of the maximum principle stresses at these three points, we find that they are related as 
C < A < B, and we expect generally that fracture arises at B before either A or C. If the 
contraction proceeds while the fixed boundary condition on the bottom is maintained, the 
lateral side breaks near the bottom before the division of the cell, and the fixed boundary 
condition can not persist. Hence we need to consider the displacement of the layer with 
respect to the bottom to deal with this problem correctly. 

In our one-dimensional model, we break a spring when its energy exceeds a critical value. 
We assume that the corresponding critical energy density is independent of both L and H. 
As mentioned above, this is equivalent to the critical stress condition in one-dimensional 
models. Representing the critical energy density with the corresponding linear shrinking 
rate Sb by the fragmentation conditions are described by the rules 




The horizontal spring is cut off; the 



(3.1) 



cell is divided. 
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I The vertical spring is cut off; the 

e 2 (x) > -hs 2 b => (3.2) 
bottom of the layer breaks. 

Here we apply the same condition to the vertical springs in order to enforce that the lateral 
side breaks before the division of a cell under the fixed boundary condition. 

We can easily determine the analytical solutions. The functional variation of the energy 
(2.4) on u(x) is obtained in the form 

5E = J dx{—k\u xx + kz(u — w)}5u 

+ [h(u x + s)Su}^ L , (3.3) 

and we obtain both the balanced equation 

u xx = q 2 {u-w) (3.4a) 

and the stress-free boundary condition 

u x + s = at x = ±L. (3.4b) 



First we assume the fixed boundary condition without slip on the bottom: w(x) = for 
\x\ < L. The solution of the equations ( |3.4j ) is then 

s sinhgx 

U{X) = - q ^L- (3 ' 5) 
The deformation almost only appears near the lateral boundaries, because of exponential 
dumping. The energy densities of horizontal e\(x) and vertical springs ei{x) are maxima at 
the center of a cell x = and at the lateral boundary x = L, respectively, and these points 
have the greatest possibility of breaking. Then energy densities are calculated as 

*<°> = 5* 14 " (' " =y (3 ' 6a) 

and 

e 2 (L) = hi lS 2 (tanhgL) 2 . (3.6b) 
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Although they both increase with shrinking, ei(0) is always less than e2 (-£/)• 

If Sb is smaller than s s , that is, if breaking occurs before slip, the breaking condition 
( p.2[) for the vertical spring on the lateral side is the first to be satisfied. To identify the 
effect of slip, we consider the fragmentation of a cell with the assumption that neither the 
slip nor the breaking of the vertical springs occurs even with the fixed boundary condition. 
For the first breaking of the horizontal springs, we apply the fragmentation condition ( |3 . 1\) 
to ( |3.6a|) and obtain the relation between the system size and the shrinking rate: 



qL = arccosh . (3.7) 

s - s b 

As indicated with the solid line in Fig. || qL drops rapidly at s/sf, ~ 1 and then vanishes 
slowly as s increases further. Because each breaking divides the system into rough halves, 
the size of a cell decreases with shrinking. The figure displays the typical development of 
the size of a cell with the stair-like function of the dot-dashed line and the arrows. Because 
the system size L is scaled by the thickness H, we see that the system is divided into a size 
smaller than H after sufficient shrinking. 

If Sb is larger than s s , the layer starts to slip from the lateral sides when e2(L) = k\s 2 s j2. 
The shrinking rate at that time is given by 

(3.8) 



tanh qL 
We next investigate this case. 

We suppose that the symmetrical slip from both lateral sides is directed toward the 
center and only consider the half region x > 0. The function w(x) becomes finite in the slip 
region x s < x < L and remains zero elsewhere, where we introduce x s as the starting point 
of the slip region. The slip displacement w(x) = w (x) is expressed by the displacement 

u(x) = uo(x) as 

< x < x s 

w (x) = { . (3.9) 

u {x) + ^ x s < x < L 
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Equation (|3.4|) then take the form 



q 2 u {x) 0<x<x s 
u 0xx = { , (3.10a) 

— qs s x s < x < L 



b.c: uo x + s = at x = L , (3.10b) 
and the matching conditions are 

w (x), u (x) and dU g!f are continuous at x = x s . (3.10c) 
We derive the solution uq(x) in each region and obtain 



A sinh qx < x < x s 



U {X) 



qs s [L - \x 



(3.11) 



x + B X, < x < L 



The three matching conditions give the integral constants A and B and yield the equation 
to determine x s , 

q(L -x s ) = — - - — ^ . (3.12) 

s s tanh qx s 

At x s = L, this reduces to (|3.8|) . This form can be approximated as L — x s ~ (s — s s )/qs s 
for qx s ^> 1 and as qx s ~ (s/s s — qL)^ 1 for qx s <ti 1. 

We calculate the energy density ei(0) again and substitute this into the breaking condi- 
tion Q3.1|) . This gives the equation 

1 > 8 X (3.13) 



s s sinh qx s s 

Eliminating x s from the equations ( p,12|) and ( p,13| ), we obtain the following relation 



between the system size and the shrinking rate at the first breaking: 



S s \ S / ( S Sfe N 2 



qL = arcsinh — + Wl + . (3.14) 

\s-s b J s s V V s s J 

We see that qL is a decreasing function of s. It decreases slowly to the limiting value Sb/s s 
after the rapid drop in the range s?, < s < s b + s s . 
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Figure 01 exhibits two curves of the shrinking rates at the start of slip (|3.8|) and at the 



first breaking ( |3.14j ), where half of the system size L is represented on the vertical axis as 
in Fig. [| After full contraction, the final size of a cell is close to the asymptotic value of 
the curve defined by ( |3.14| ), qL ~ s&/s s . The region without slip also becomes smaller, and 
its final size is given by qx s ~ s s /(s — Sf,), where we assume qx s <C 1 in fl3.12| ). With the 
original scale, we obtain 

L ~ and x s ~ — for s — Sb ^ s s . (3.15) 

q s s q s-s b 

Thus L and x s are proportional to the thickness of a layer H, although there is the possibility 
for them to be modified through s s if the frictional force depends on H. 



The first equation of (|3.15|) is consistent with the experimental results of Groisman 
and Kaplan [fj] for the final size of a cell after full desiccation, as mentioned in Sec. |. 
The assumptions used in this analysis are also consistent with those in their qualitative 
explanation, where they considered the balance between the frictional force and elastic force 
0. In addition, when we peel the layer of an actual mixture after drying, we often observe 
a circular mark at the center of each crack cell on the bottom of the container. Its size is 
approximately equal to the thickness of the layer. We can understand these marks as the 
sticky region |x| < x s . 

B. The Griffith Criterion 



Next we apply the Griffith criterion [£9|] to the entire system as the fragmentation con- 
dition in the place of the critical stress condition. This was used by Sasa and Komatsu in a 
different model 0. 

First we again assume the fixed boundary condition, where neither the slip nor the 
breaking of the vertical springs occurs. The Griffith criterion introduces the creation energy 
of a crack surface per unit area V and assumes the cracking condition that the sum of the 
creation energy and the elastic energy decreases due to breaking. We write the elastic energy 
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of a system — L < x < L as E(2L). We consider the case in which the cell with size 2L (the 
system size) is divided into exact halves. The alternative fragmentation condition to ( |3 . 1| ) 
is given by the equation 

AE(2L) = E(2L) - 2E(L) > TH. (3.16) 

We calculate ( f2.4| ) by using ( |3.5| ) to obtain the elastic energy E{2L) for the fixed boundary 
condition. With the original scale, it is given by 

E(L) = k lS 2 HL ( 1 - — tanh ^ ) . (3.17) 
\ qL H J 

As a result, we obtain the following relation in the place of (|3.14 ) for the shrinking rate 
at the first breaking: 



Here, 



A£(L)=(f) 2 and s T = X^. (3.18) 



AE(L) S ff^ = 2 tanh ^ - tanh ?± 



1 gL > H 

(3.19) 

J (f ) ^ « H 

The corresponding curve is indicated with the dotted line in Fig. where the shrinking rate 
s is scaled by Sp. This curve agrees quite well with the solid line representing the previous 
results ( [3.7|) , so we again find that cells are divided into a size smaller than H after breaking. 
We however note that sp depends on the thickness H, although both and sr represent the 
shrinking rate at the first breaking for an infinite system. Because the ratio of the surface 
energy T to the elastic constant is a microscopic length for ordinary materials, sr is inferred 
to be very small. Therefore, with the Griffith criterion, we usually expect that sr is smaller 
than s s and no slip occurs before breaking. 



Next we show that, even if sr is larger than s s , the Griffith criterion does not yield the 
proportionality relation of the final size of a cell to the thickness of the layer H. Elastic 
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energy is consumed not only by the creation of the crack surface but also by the friction due 
to slip on the bottom. We again consider the breaking of the system (— L < x < L) into 
exact halves. An alternative Griffith condition is given by 

AE S (L) = E S (2L) - 2[E' S (L) + W s ] > TH, (3.20) 

where E S (2L) and 2E' S (L) represent the elastic energies of the system before and after 
breaking, respectively, and 2W S is the work performed by the frictional force due to slip. 

As the state just before breaking, we consider a cell with symmetric slip regions. This 
state has been derived in ( |3.9|) , (|3.11|) and ( [3.12| ). The elastic energy E S (2L) is obtained by 



calculating ( |2.4|) in the form 

k\S 2 s 



E S (2L) 



Q 



\q\L - x s f + (-) qx s - — 
3 \s s J s s 



(3.21) 



where the width of the slip region L — x s is determined as a function of L and s/s s by fl3.12| ) 



In order to estimate W s and E' S (L), we need to investigate the detailed process of frag- 
mentation. Here we imitate an actual quasi-static fracture in two dimensions by using a 
hypothetical quasi-static process in the one-dimensional model. We introduce a traction 
force on crack surfaces which prevents the crack from opening and obtain the final state of 
this process with the stress-free boundaries by stipulating that the strength vanishes quasi- 
statically. The work of the hypothetical traction force is considered to be the opposite of 
the creation energy of the crack. Let us imagine the right half < x < L just after the 
breaking at the center x = 0, where the traction force works at x = to the left. Because 
of the relaxation of the traction, the slip region x s < x < L before the breaking vanishes 
immediately. As the traction decreases, a new slip region is created in < x < x r on the 
side of the crack. If the contraction ratio s is much larger than s s , we may assume that x s is 
smaller than x r at the end of this process, because qx s < 1, and the new slip region expands 
to x r ~ L/2. The slip displacement w(x) in the state is given by the initial condition (|3.9|) 
and the slip condition (|2.5|) as 
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The solution of 



W(X) 



IS 



Wq(x) X r < X < L 

u(x) — — < X < x r 



(3.22) 



u(x) 



Uq(x) + C[ cosh q(x — L) x r < x < L 



~^qs s x sx ~\~ C2 



< X < X r 



where the conditions of the continuity of w(x), u(x) and u x (x) at x 
constants C[ and C 2 and produce the equation for x s : 

q(L — 2x r ) = 2tanhg(L — x r ). 



(3.23) 



x r determine the 



(3.24) 



We calculate the elastic energy ( |2.4j ) from ( |3.22| ) and ( |3.23| ) and obtain the energy at the 
end of this process, 



2K(L) = ^ {\<?[*> + (L - x r f] - qL) . 



(3.25) 



Because slip occurs with a constant frictional force from the previous assumption, the work 
W s can be expressed by the integral of the total distance of slip, 

F s r L 



a Jo 



dx\w(x) — Wq(x)\, 



(3.26) 



and it is calculated from ( |2.6|) , (p.22|) and ( p.23| ) as 



2We 



2 r 



-q 6 {x 6 s - 2x 6 r ) + q"Lx 2 r - IqL - — j // 



2„2 



(3.27) 



In order to know the scaling relation of the final size of a crack cell after full desiccation, 
we assume qL ^> 1 and the limit of the full contraction: s/s s — ► 00. Because Q3-12 ) 
and ( p.24 ) give the approximate equations qx s ~ s s /s 1 and x r ~ L/2, respectively, 
(CT),fl£2TD, fl£25D and (|3~27| ) result in the equation 



(3.28) 



With the original scaling, the Griffith criterion ( 3.2(Jp gives the scaling relation of the final 
size of a cell for the thickness H, 
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2 

qL > ^6H [j-j 3 oc Hi, (3.29) 



where we use sr defined in ( |3.18j ), and the condition sr ^> s s is necessary from the assumption 
qL 1. Thus the Griffith criterion gives the different scaling relation because of the 
dependence of sr on H, although we obtained the proportionality relation ( |3.15| ) under the 
critical stress condition. 



As a result, the critical stress condition and the Griffith condition lead to different 
relations between the final size of a crack cell and the thickness of a layer. Experimental 
results seem to support the former. Although the results should be discussed further, of 
course, they suggest the possibility that the dissipation in the bulk can not be neglected for 
the fracture of a mixture of granules and water. 



IV. THE DEVELOPMENT OF CRACKS IN TWO DIMENSIONS 

The one- dimensional model we have discussed to this point idealizes the process of crack- 
ing, treating cracks as one- dimensional structures forming one at a time parallel to one an- 
other. In order to consider the development of cracks and their pattern formation, we must 
extend this model to two dimensions and include the relaxation process of the elastic field. 

Although mixtures containing granular materials that are rich with water generally pos- 
sess viscoelasticity, visible fluidity can not be observed at the time of cracking after the 
evaporation of water with desiccation. Therefore, we assume that only the relaxation of 
strain contributes to the dissipation process in a linear elastic material. For simplicity, we 
assume that the system has the only one characteristic relaxation time. 



In the one-dimensional model, the total energy (7a) consists of terms representing a 
one-dimensional linear elastic material and the vertical shearing strain, the both of which 
are quadratic in the displacements. Expanding the elastic material to the horizontal xy 
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plain, we naturally obtain the extended energy in two dimensions as 

E = J dxdy{e x (x, y) + e 2 {x, y)}, (4.1a) 

1 / 1 \ 2 

ei(x, y) = -n(u u + C) 2 + /x [u ik - -uu6 ik J (4.1b) 

and e 2 {x,y) = -k 2 (u - w) 2 , (4.1c) 

where Uij = (uij + Uj t i)/2 and Uij = diii/dxj. In analogy to the one- dimensional model, 
the two-dimensional vector fields u(x, y) and w(x, y) represent the displacement of a layer 
from the initial position and the slip displacement on a bottom, respectively. The space 
coordinates x and y and the displacements u(x,y) and w(x,y) are again scaled by the 
thickness of a layer, H. The expression in ( |4.1b| ) is the energy density of the two-dimensional 
linear elastic material with a uniform free surface-shrinking rate C. 

The shrinking speed C can be neglected from the assumption of quasi-static contraction. 
The time derivative of (|4.1|) is given by 



E = J dxdy{[-a ijd + k 2 {ui - w t )] 



-Wui-wW + fdSnwii, (4.2a) 
<7ij = (Xuu + nC)8ij + 2/iUij and k = A + fi, (4.2b) 

where § dS represents the integral along the boundary of the cell and n is its normal vector. 

The dissipation of energy arises from the non- vanishing relative velocities of neighboring 
elements in a material, and then the time derivative E can also be represented as a function 
of them. As is well known, E can be written in a form similar to the energy due to certain 
symmetries [^UJ. We have 

2 



E = - J dxdy{K'u 2 u + 2// (ii ik - -uuS ik 



+K(u-w) 2 } 
dxdy{[-a' ij:j + k' 2 (ui - Wi)\ui 

—k' 2 (ui — Wi)wi} — j> dSuju'^Ui (4.3a) 
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to the second order, where 

(jjj = X'uiAj + 2ii'u i:j and k' = X' + \i . (4.3b) 

Although the constants k', /jl' and k' 2 are generally independent of k, /z and k 2 , we assume 
they take simple forms with one relaxation time r, writing a'^ = rdaij/dt, or equivalently, 

k' = tk, /jf = r/i and k' 2 = r/c 2 . (4.4) 

Equations (f|7|), (|4.3j ) and ([O]) yield the time evolution equation of u(x,y), 

(l + - fc 2 ( Wi - Wi)} = 0, (4.5) 

and the free boundary condition, 

(l + r^j a tJ n 3 = 0. (4.6) 

As mentioned above, the shrinking rate C only appears in the free boundary condition. 

Except for the effect of shrinking and slip, the above equations are essentially the same 
as the Kelvin model, proposed for viscoelastic solids. Because the existence of a bottom 
causes a screening effect through the term k 2 Ui, the elastic field decays exponentially in the 
range of the thickness of a layer, i.e, the unit length in (|4.5|) . Here the stress in the material 
is (1 + Td/dt)a,ij by adding the dissipative force. With the definitions 

£/i= (l+rJ^Ui, W^^l + T^jw, (4.7a) 

and 



% = (l + ^K> ( 4 - 7b ) 



5|) and (14.61) can be rewritten as 



E yJ = k 2 {Ui - Wi), (4.8a) 
= (XU U + KC)8 i5 + 2nU iS , (4.8b) 
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doing with the free boundary condition 

UjZij = 0. (4.8c) 

Therefore, [/, satisfies the balanced equations of an ordinary elastic material without dissi- 
pation. 

We expect that the state of the water bonds in a mixture can be represented by cr^, 
i.e. the stress excluding the dissipative force, rather than by the stress Ejj itself because 
Oij is a function of the strain Uij. For this reason we introduce the breaking condition by 
using Oij in the following analysis. The propagation of a crack in the Kelvin model has been 
investigated by many peoples pT]- p5| , p^] . Although the stress field diverges at a crack tip 
in the continuous Kelvin model, it is possible that the divergence of is suppressed by 
the advance of a crack. We calculate the elastic field around a crack tip for a straight crack 
which propagates stationarily in an infinite system. Because near the tip of a propagating 
crack there is little slip, as the simulations in the next section indicate, we can assume 
w ~ 0, and then Wi(x, y) = in ( }4.8| ). Although our model is incomplete in the sense 
that the divergence of the stress can not be removed in continuous models of a linear elastic 
material, we expect that the following discussions are valid. 

We consider a straight crack with velocity v that coincides with the semi-infinite part of 
the x-axis satisfying x < vt in two-dimensional plane. The stress field satisfies the stress-free 
boundary conditions on the crack surface, and the displacement u vanishes as \y\ — > oo. We 
define the moving coordinates (£, y) as £ = x — vt and assume reflection symmetry on the 
x-axis. The boundary conditions are given by 

12 yy = on y — 0, £ < 
U y = ont/ = 0, £>0 

(4.9) 

= on y = 
Ui — > for \y\ — > oo 
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The stress field under the above boundary conditions can be obtained by the Wiener-Hopf 
method. Fortunately, this problem reduces to the following solved problem for the mode 
I type of a crack. We consider a stationary crack along the negative x-axis in completely 
linear elastic material without contraction, where we include the inertial term with the mass 
density p. When a uniform pressure a* is added on the surface of the crack from time t = 0, 
the stress field it® is given by the equations 

d 2 u° 



a 



dt 2 



and = A<<% + 2pu% 



(4.10) 



and the boundary conditions 



< 


= -a*e(t) 


on y = 


0. 


x < 


<- 


= 


on y — 


0. 


x > 


< 


= 


on y — 









-> 


for \y\ 




OG 



(4.11) 



These become identical to the previous equations when we apply the Laplace transformations 
on time, 



Ui(x,y,r]) = / dtu° i (x,y,t)e * 
Jo 



(4.12a) 



and 



^i(x,y,T})= dto?Ax,y,t)fe~ 



■rjt 



(4.12b) 



and make the replacements rj 2 p = a* = kCtj, = £° + nC5ij and a; — > £. Here rj can 
be taken equal to 1 because the correspondence holds for any rj. 

Using the analytical solutions [|32] of ( |4. 10| ) and ( |4.11| ), we shall consider the stress in 
front of the crack, i.e. £ yy (£, 0) where £ > 0. The above replacements give the solution of 
our problem as, 

j roo—Oi 



0) = kC — / rfC Jm£°(-C) e~ CC + 1 



a: 



and E0(C)4f§#-l^ 



(4.13) 
(4.14) 
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where F + (() is an analytical function in the region Re( > —a such that 



F+(C)-r s for |C| ^ oo 



(4.15) 



and 



*+(0) 



2/z(A + ji) 



(4.16) 



\| a(A + 2//) 2 ' 

Because 7mS°(- Q approaches — F + (0)/\/C as ICI increases, it can be approximated in the 
region C 1 around the crack tip by the equation 



e w (£, o) ~ kC f dcch-x + lj 



We see that the stress T, yy diverges in inverse proportion to the square root of the distance 
near the tip. 

Solving (|4.7b ) in the moving system £ = x — v t, we get 



1 r°° £-£' 
a yy^ = —j i d^ yy ^',0)e— , for£>0. 



(4.18) 



and the following equation is obtained by substituting ( |4.13| ) into this equation: 

1 roo-w I m £°(-() _ c 



<7yy(Z,Q) = *C S / d C r \ < 

I 7T Ja-Oi rt»C + 1 



-a 



n. 



(4.19) 



We first examine the behavior of a yy in the region a£ ^> 1 far away from the tip. Because 
the integral in Q4.19D can be approximated as 

^ 7mS°(-0 e _ ce 

a-Oi TV( + 1 



o tv(Q + a) + 1 

Imt°(-a) e" Q? 
rva + 1 £ ' 

cr„„ decays exponentially to the uniform tension as 



(4.20) 



<7 W (£, 0) ~ kC (j e '^ + A for «C > 1, 



(4.21a) 
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where 

D = _ Im ^(-a) 
n(rva + 1) 

The characteristic length of the decay is H/a with the original scale, which is of the same 
order as the thickness H for an ordinary elastic material. 

In the region C 1 near the tip, we approximate the integral with the asymptotic form 
of JtoE°(— Q for large ( as 

1 r ~° V /mS°(-Q c _ C£ 



poo—w 
" - / d( 

71 Ja-Oi TVQ + 1 

POO 

~ / d( — -c 

Jo 



7t(tv( + 

^Be^erfc [JL I . (1.22a) 

\J 7TTV \ V TV 



roo 2 ! 2 

erfc(x) = / dte^ ~ — e~ x for x — > oo, (4.22b) 

ji 2x 



7T 



and erfc(O) = ^— , (4.22c) 



and then we obtain the equations 



<t w (£,0)~{ / fora£«l. (4.23) 



With the original scale again, <7ij is almost proportional to y H/£ in the middle region from 
tv to H/a around the crack tip, as is the case for the stress E i3 -. However, is bounded 



at the tip by the proportional value to y 1 ' H/rv. 



Thus the stress excluding the dissipative force, i.e. a^, is kept from diverging by the 
movement of the crack tip. Therefore we can introduce a critical value a th as a material 
parameter again and assume the breaking condition at the tip 

lkn <7 w (f , 0) > a th = KC th , (4.24) 
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where the constant Cth represents the shrinking rate corresponding to the critical value. 
For a stationary propagating crack, we find the velocity v by substituting cr yy (^, 0) into the 
above equation in (|4.23|) as 




(4.25) 



Although this equation is not valid near the sound velocity because we have neglected inertia, 
we expect that a material cracks at a shrinking rate C below Cth due to inhomogeneity. Thus 
( |4.25| ) explains why the propagation velocity observed in experiments is very small compared 
to the sound velocity. In addition, it suggests the proportionality relation between the 
thickness and the propagating speed, which should be experimentally observable. 



Next we note the validity of Q4.25 ) for very slow speeds. Although the divergence of a 



mi 



is suppressed by the advance of a crack, as we see in ( |4.23|) , the size of the screening region 
is approximately rv. Because particles of size R = 0.01 ~ 1mm are used in experiments, the 
above breaking condition fl4.24p is available for velocities v 3> R/t, where the continuous 
approximation is valid. For very slow velocities, rv <C R, for example, it may be possible 
for the defects in a material to arrest the growing of cracks. 

The speed of a real crack measures v ^ 2mm /min for the thickness of a layer of coffee 
powder, H ~ 6mm, as observed in experiments ||. Although we have not yet specified 
the origin of the relaxation, we attempt to estimate the relaxation time r arising from the 
viscosity of the water in the bonds among particles. Supposing that the diameter of a particle 
of the coffee powder is about R ~ 0.5mm and the viscosity of the water is v ~ 1mm 2 /s, we 
obtain r ~ R 2 jv ~ 0.25s and tv/R ~ 0.02 < 1. This rough estimate suggests the possibility 
that the above continuous approximation is imperfect in the case of a relatively thin layer. 

Groisman and Kaplan reported the interesting experimental results that the propagation 
speed exhibits a wide dispersion even among cracks growing at the same time, although the 
speed of each individually is almost constant on time. It is a future problem to understand 
the relation of the dispersion to the inhomogeneity of materials. 
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V. THE PATTERNS OF CRACKS 



Cracks appear one after another with shrinking, and spread over the system to create 
a two-dimensional pattern. In this section, we report on a study of the formation of crack 
patterns using numerical simulations of the two-dimensional model introduced in the pre- 
vious section. We make the natural extension of both the breaking condition and the slip 
condition employed in the one-dimensional model by using the energy densities defined in 
the microscopic cells of the lattice in the simulations. In the two-dimensional model, the for- 
mer is simpler than the critical stress condition because it neglects the direction of both the 
stress and the microscopic crack surfaces. In addition, the slip condition implicitly assumes 
a sufficiently short period of slip in comparison to the relaxation time of the elastic field 
because of the balance between the frictional force and elastic force. We report the results 
of our simulations after describing the discrete method and these extended conditions. 



As is the case with most fingering patterns, the growing of cracks is influenced strongly 



by the anisotropy of the system. We used random lattices m our simulations to 

consider uniform and isotropic systems in a statistical sense. 

Many fracture models employ a network of springs or elastic beams to model an elastic 
material |3^- ^5[ , |I1] , |I3] ,^|, |ID| , ^5[ , ^D| . However, it is generally difficult to calculate the elastic 
constants for an elastic material modeled by such a lattices. Therefore, instead of such 
networks, we consider each triangular cell in a random lattice as a tile of the elastic material 
with uniform deformation. A fracture is realized by removing any cell whose energy density 
exceeds a critical value, as we explain below. 

We construct a two-dimensional model as is illustrated in Fig. [5[ Each site of the lattice 
is connected to an element on the bottom with a vertical spring similar as that used in the 
one-dimensional model. Figure ^| displays a part of the random lattice which represents a 
horizontal elastic plane. The random lattice is composed of random points to form a Voronoi 
division of them. We number the sites and the triangular cells in the lattice and express 
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the area of the Voronoi cell around the nth site as V n (n = 1, 2, 3, N) and the area of the 
mth triangle as T m (m = 1, 2, 3, Nt)- 

The displacement of the nth site u^(t) is obtained from the U^(t) by the definition 
( |4.7a| ). A simple Euler method gives the following equation with discrete time At: 

U W(t + At) = u^ n \t) + — (U (n) (t) - u< n >(*)). (5.1) 

r 

Because U(£) satisfies the balanced equation of an ordinary elastic material at any time, 
it minimizes the energy E defined by (O) through the replacements u — > U and w — > W. 
E consists of the energies of both the vertical springs and the horizontal elastic plain. We 
calculate the latter by summation of the energies of the triangular cells in the random lattice, 
where the mth triangular cell is assumed to consist of a linear elastic material with uniform 
strain tensor . Thus we obtain the equations 

E = Y^T m e i ? ) + Y.V n ef\ (5.2a) 

m=l n=l 

e! m) ee I«(17<™> + Cf + u (u^ - \li m) S S k) a (5.2b) 
and 4 n) = ^ 2 (U (n) - W^) 2 , (5.2c) 

where and W^ n ^ represent the displacement of the nth site and the slip displacement 
of the nth element along the bottom, respectively. Although the rule of repeated indices is 
applied to j, k and I, as usual, the summations over m and n are expressed by the symbol 
J2, where J2' m represents a summation that excludes broken triangle cells. 

The quantity is the elastic energy of the mth triangle cell which is calculated from 
u\™\ For the following explanation, we express the vertices of the mth triangle as n = 1, 2, 3 
and their initial equilibrium positions as x< n > = (x (n) ,y (n) ), as is shown in Fig. |5|. Assuming 
uniform deformation in the triangle, the strain tensor U^f 1 = \{Uy^ + Uj^) ls given by the 
displacements of the vertices U^) as the equations 

U£ ] = ^HVMUW, (5.3) 
U™ = ~e ijk x^U^, (5.4) 
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= W *kh™uW - xMuPl (5.5) 

where 

T = ei jk x®yM, T m =~\T\, x^=x«-x^, (5.6) 

and e.ijk is Eddington's e. These equations are easily obtained from the first order Taylor 
expansions of U (n) = U(xW,t) in the triangle. Thus we can calculate the energy (|5.2|) on 
the random lattice and obtain from its minimum. 



A fracture is represented by the removal of triangle cells, not by the breaking of bonds, 
in this model. This gives a direct extension from the one- dimensional model, although 
it neglects the microscopic direction of the stress and the cracking in a triangle cell. We 
calculate the elastic energy density of the mth triangular cell from the true displacements 

, and assume a critical value for the breaking, which is represented by the corresponding 
shrinking rate Cf, as Thus the breaking condition is given by 

, s i The mth triangle cell is re- 

e! m) > -kCI (5.7) 
A moved, 

where 

e! m) = \*(uP + C) 2 + u (v$ - ^4 m) ^) 2 , (5-8) 

and u& is calculated from using the method explained above for Uj™\ 

The slip condition for the elements on the bottom is also similar to that in the one- 
dimensional model. We introduce the maximum frictional force as a constant and assume 
the balance of the fictional force against the sum of the elastic force and the dissipative 
force, i.e., \k 2 (1 + rd/dt) (u — w)| = k 2 \XJ — W|. Therefore, the slip condition for the nth 
element can be written by using as 

W( n ) is moved along the 
e>2^ > ^ f° rce t° the position where (5.9) 
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The slip displacement w( n ) is calculated from W*™) by the definition ( |4.7a|) as 

wW(t + At) = wW(t) + — (W (n) (t) - w (n) (t)). (5.10) 

r 



We carried out the numerical simulations of the model by increasing the shrinking rate 
C in proportion to time t with a constant rate C. We repeated the following procedures at 
each time step: 

1. The contraction rate increases as C(t+ At) — C(t) + AtC. 

2. The displacement u is calculated from U which is the function minimizing the energy 

D- 

3. The slip w is calculated from W, which is given by the conditions (|5.9| ) at the all sites. 

4. If some triangular cell satisfy the breaking condition ( |5.7|) , it is removed. Its energy is 
not included in subsequent calculations. 

If more than one cell satisfies the breaking condition at step 4, we repeat the step 1-3 using 
a smaller time step. 

In the above calculations, we can take the parameters k, k 2 , Cf, and C to be 1 with loss 
of generality by scaling space, time, energy and the shrinking rate as follows, 



1 -y2 p . . £V 



E-^-kC£E } t^^t, and C = C b C. (5.11) 
2 C 

Here we write the scaled shrinking rate as C. As a result, three independent parameters 
remain explicitly in the equations: 

/*=-, t = —t, and C s = -^-. (5.12) 

Although the properties of the lattice influence the results, we used an identical lattice 
in all of our simulations, except the last in which we consider the effect of a random lattice. 
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To prepare the sites in the random lattice, we arranged points in a triangular lattice with 
mesh size 0.01 inside a square region lxl and shifted the x, y coordinates of each point 
by adding uniform random numbers within the range ±0.005. Because their distribution is 
almost uniform and random, we connect them by Voronoi division to make the network of 
the random lattice. Then the square region is extended to the size L x L, which is related 
to the original system size L as 

L.M (5,3) 



We use the conjugate gradient method El] with a tolerance 10 6 to find the minimum 



points of functions on free boundary conditions. The time step At is changed automatically 
in the range to an upper bound 10~ 4 — 10~ 3 . In our simulations, at most one cell is removed 
in any given timestep. In contrast, we note that in the model without a relaxation process, 
a crack propagates instantaneously at a fixed shrinking rate, because the breaking of a cell 
necessarily changes the balanced state of the elastic field and often causes a chain consisting 
of the breaking of many cells occuring simultaneously. 

We show the typical development of a crack pattern in our model in Fig. [7], where the 
parameters are L = 10y/2, ft = 1, f = 0.01 and C s = 0.5, and the black area indicates 



the removed triangular cells. The gray scale represents the energy densities (|5.8|) in Figs. 
[7| (a),(b) and (d), and each dot in Fig. [7| (c) is a slipping element under the condition 
( |5.9| ). (7=1 (i.e. C = C^) corresponds to the shrinking rate at the first breaking for an 
infinite system. The first breaking in the simulations occurs slightly below C = 1 for most 
parameters because of the randomness in the lattice. At that time, no slip has yet begun in 
the most of the system, except near the boundaries. 

After this, some crack tips grow simultaneously in the whole system, and the crack 
pattern almost becomes complete at the shrinking rate near C = 1. Here we see the white 
circular marks around the center of the crack cells, as shown in Fig. 0-(c). These represent 
the sticky regions without slip. Similar marks can be observed on the bottom of a container 
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in actual experiments, as we mentioned in Sec. III. 



Figure |8] graphs the development of the total energy ( |4.1| ) with shrinking for the three 
cases f = 0.1, 0.01 and 0.001. The energy increases with contraction. However, it is released 
and dissipates due to successive breaking and becomes almost constant with increasing 
shrinking rate. Our simulations were carried out until C = 10. The crack patterns changed 
little when C becomes large, while the circular marks shrank gradually. 

For fast relaxation, new cracks grow from the lateral side of another crack almost per- 
pendicularly. Figure || displays a snapshot of a crack pattern for the very small relaxation 
time f = 0.001. As r becomes smaller, the cracks tend to propagate faster and grow one 



at a time, in agreement with the assumption in the one- dimensional model. In Fig. [H| 
we compare how the total number of broken triangular cells increases with shrinking for 
f = 0.1, 0.01 and 0.001. The number of broken cells is almost proportional to the total 
length of cracks. It increases like a step function for f = 0.001. This indicates that the 
cracks are formed one by one. 



For slow relaxation, in contrast, the growing cracks from fingering-type patterns [3"5-3"5 
such as similar to those seem in viscous fingering. As is shown in Fig. many cracks tend 
to grow simultaneously from the center toward the boundaries. They are accompanied by a 
series of tip splittings and the total length of the cracks increases smoothly with contraction. 
As r becomes larger, it takes more time to complete the crack patterns because of the slower 
propagation of cracks. 



We can see the influence of slip on the crack patterns by changing C s with the other 



parameters fixed. Figures 12 and O are snapshots of a crack pattern after full shrinking: 



C = 10 for C s — 1.0 and C s = 0.1, respectively. Figure |T3| shows the total number of broken 
triangular cells for the various values of C s . This figure indicates that the crack patterns at 
C = 10 are close to final states. As we expect, the final size of a cell becomes larger with 
smaller C s . Figure [15| plots the total numbers of broken triangles at C = 10 for C s . They 
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increase monotonously in this range with an almost constant rate. 



Next we change the elastic property with fx. From the equations of <Jij ( |4.23| ) and the 
crack speed ( |4.25| ) in Sec. [TV], we expect that, as fx becomes smaller, the stress, excluding 
the dissipative force, has a weaker concentration at a crack tip, and the cracking speed is 



smaller. Figure [16] shows a snapshot of a crack pattern for fx = 0.02. We find that the cracks 
become irregular and jagged lines and that they propagate slowly. The crack patterns also 
reach completion more and more slowly as fi decreases, as is shown in Fig. [17|. 

All of the above simulations were executed on the same random lattice. For comparison, 
we also performed a simulation using a regular triangular lattice in the place of the random 
lattice. Figure [18] shows the result using the same parameters used in Fig. ^ except for the 
difference of the lattices. It is clear that the anisotropy of the triangular lattice is reflected 
in the direction of cracks. 



Thus we can reproduce patterns similar to those of actual cracks using our two- 
dimensional model. The dependence of the formation of cracks on the relaxation time 
and the elastic constants should be compared with actual experiments in more detail. The 
experimental results of Groisman and Kaplan suggest a transition in the qualitative nature 
of patterns as the thickness is changed. This is point (3) mentioned in Sec. |. Similar 
changes of patterns are observed for slow relaxation or small fx in our simulations. However, 
our model does not contain the thickness H explicitly because of the scaling with H, and we 
have no experimental data for the dependence of the other parameters on H. In addition, 
it is possible that the inhomogeneity in a material plays an important role in this change 
because the size of a particle can become significantly large if H is made sufficiently small 
0. Obtaining a more detailed understanding that address these points is left as a future 
project. 
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VI. CONCLUSIONS 



We studied the pattern formation of cracks induced by slow desiccation in a thin layer. 
Assuming quasi-static and uniform contraction in the layer, we constructed a simple model 
in Sec. O. It models the layer as a linear elastic plane connected to elements on the bottom 
and considers the slip with a constant frictional force. 

In Sec. [IT], we considered the critical stress condition by introducing a critical value of 
the energy density. This model explains the proportionality relation between the final size 
of a crack cell and the thickness of a layer and the experimental observations on the effect 
of slip. We also considered the Griffith criterion as an alternative fragmentation condition. 
However it predicts qualitatively different results for the final size of a crack. Although these 
results need more detailed considerations, they suggest the possibility that dissipation in the 
bulk may be important in these materials. 

In Sec. [iy|, we extended the model to two dimensions and introduced the relaxation 
of the elastic field to describe the development of cracks. This is essentially the same as 
the Kelvin model for viscoelastic materials. Because the stress excluding the dissipative 
force does not diverge at the tip of a propagating crack, we introduced a critical value as 
the breaking condition in front of a moving crack. Assuming the existence of a stationary 
propagating crack, we obtained an estimation for the crack speed in closed form within a 
continuous theory. This estimation explains the very slow propagation of actual cracks and 
predicts the proportionality relation between the crack speed and the thickness of a layer. 
It is an open problem to understand the origin of the dissipation we introduced intuitively 
and the role of inhomogeneity in the stability of a crack. 

In Sec. 0, we carried out numerical simulations of the two-dimensional model to investi- 
gate the formation of crack patterns. By using the energy density defined on the microscopic 
cells of a lattice, we introduced a simplified breaking condition from the direct extension of 
the one- dimensional model. We used a random lattice to remove the anisotropy of the lattice 
and obtained patterns similar to those observed in experiments. We found that cracks grow 
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in qualitative different ways depending on the ratio of the elastic constants and the relaxation 
time. It is important in the connection with fingering patterns that for the slow relaxation 
crack patterns are formed by a succession of tip-splittings rather than by side-branching. 
We need a better understanding of the role of inhomogeneity to explain the transition in 
the nature of patterns as a function of the thickness reported in the experiments. 

For the experimental results of Groisman and Kaplan, which we mentioned in the points 
(l)-(3) in Sec. |, we believe the present results give qualitative explanations for (1) and (2) 
and a clue for (3), although more considerations is necessary. 
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FIGURES 
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FIG. 1. The one-dimensional model 
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FIG. 3. The shrinking rate at the first 
breaking under the fixed boundary condi- 
tion is plotted function of the system 
size, where the vertical axis represents qL 
scaled by H and the horizontal is s scaled 
by Sb or sr ■ The dotted line represents the 
result for the Griffith criterion. 



iB 



FIG. 2. The cross section of a shrink- 
ing cell of an elastic material. 
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FIG. 4. The shrinking rates at the 
start of slip and at the first breaking are 
plotted for the system size, as in Fig. H. 




FIG. 5. The two-dimensional 
model. The blank triangles repre- 
sent a crack. 




FIG. 6. A part of a random lattice 
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(c) C = 1.46 (d) C = 10.00 

FIG. 7. The time series of a crack pattern for L = 10\/2, ft = 1, f = 0.01 and C s = 0.5. The 
black area represents cracks. The gray scale in (a),(b) and (d) indicates the energy densities of the 
triangular cells of a lattice, and the dots in (c) indicate slipping elements. 
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1 2 
contraction ratio C/C b 

FIG. 8. The change of the total energy 
with contraction for L = 10\/2, fi = 1, 
C s = 0.5 and f = 0.001 for the solid line, 
f = 0.01 for the dotted line, and f = 0.1 
for the dot-dashed line. 
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FIG. 10. The change of the total num- 
ber of broken triangular cells for the sim- 
ulations of Fig. 




FIG. 9. The time development of 
cracks with fast relaxation. L = 10\/2, 
jX = 1, C s = 0.5, f = 0.001 and C = 1.54. 




FIG. 11. The time development of 
cracks with slow relaxation. L = 10\/2, 
f t = l i C s = 0.5, f = 0.2 and C = 1.09. 
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FIG. 12. A final crack pattern on a 
sticky bottom. L = 20\/2, fj, = 1.0, 
f = 0.01, C s = 1.0 and C = 10. 
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FIG. 14. The change of the total 
number of broken triangular cells for 
L = 20v / 2, p, = 1.0, f = 0.01 and 
C s = 0.1, 0.2, 0.5, 0.7, 1.0. 



FIG. 13. A final crack pattern on a 
slippery bottom. L = 20\/2 ; A = 1-0, 
f = 0.01, C s = 0.1 and C = 10. 
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FIG. 15. The final number of broken 
triangles, where we plot the values at 



C = 10 in Fig. 14 for C t 
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FIG. 16. A crack pattern for small £i. 
L =10V2, A = 0.02, C s = 0.5, f = 0.01 
and C = 10. 




FIG. 18. A crack pattern on a regular 
triangular lattice for the same parameters 
as in Fig. [7| (d). 
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FIG. 17. The change of the total num- 
ber of broken triang les for L = 10^/2, 
C s = 0.5, f = 0.01 and fi = 0.02, 0.1, 
1.0. 



42 



